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C. Weeks and M. Franz 
Department of Physics and Astronomy, University of British Columbia, Vancouver, BC, Canada V6T IZl 

We explore the possible particle-hole instabilities that can arise in a system of massless Dirac 
fermions on both the honeycomb and 7r-flux square lattices with short range interactions. Through 
analytical and numerical studies we show that these instabilities can result in a number of interesting 
phases. In addition to the previously identified charge and spin density wave phases and the exotic 
'quantum anomalous Hall' (Haldane) phase, we establish the existence of the dimerized 'Kekule' 
phase over a significant portion of the phase diagram and discuss the possibility of its spinful 
counterpart, the 'spin Kekule' phase. On the vr-flux square lattice we also find various stripe phases, 
which do not occur on the honeycomb lattice. The Kekule phase is described by a Z3 order parameter 
whose singly quantized vortices carry fractional charge ±e/2. On the 7r-flux lattice the analogous 
dimerized phase is described by a Z4 order parameter. We perform a fully self-consistent calculation 
of the vortex structure inside the dimerized phase and find that close to the core the vortex resembles 
a familiar superconducting U(l) vortex, but at longer length scales a clear Z4 structure emerges with 
domain walls along the lattice diagonals. 



I. INTRODUCTION 

Electrons hopping on the two-dimensional honeycomb 
lattice and the square lattice with a half magnetic flux 
quantum piercing each plaquette (7r-flux lattice) exhibit, 
near half filling, a linearly dispersing excitation spectrum 
characteristic of massless Dirac fermions. The honey- 
comb lattice is realized in single layer graphen^ while the 
TT-flux square lattice has the potential to be realized in 
artificially engineered semiconductor heterostructures.'^' 
Such massless Dirac fermions exhibit a host of fascinat- 
ing properties which underlie much of the current interest 
in graphene and the related systems. 

Although Dirac fermions in graphene are intrinsically 
massless, it is interesting to contemplate the effects on 
its electronic properties of a bandgap that would give 
rise to massive Dirac fermions. Experimentally, such 
a bandgap can be realized by placing graphene on a 
specific substratd^ that breaks the sublattice symmetry. 
Even more interesting is the possibility of the interaction- 
driven gap, resulting in a Mott insulating behavior. Aside 
from the conventional charge density wave (CDW) and 
spin density wave (SDW) instabilities, the structure of 
graphene allows for more interesting phases such as the 
quantum anomalous Hall (QAH) phase discussed by 
Haldane* and the quantum spin Hall (QSH) phase .EEI 
These ph ases are characterized by nontrivial topological 
invariantJ^Sini and exhibit the quantum Hall effect in ab- 
sence of magnetic field and the quantum spin-Hall effect, 
respectively. Although, as far as we know, these phases 
do not occur in natural graphene, the QSH phase was pre- 
dicted to occui^ and subsequently observed in HgTe 
quantum wells. Also, theoretical studies of these phases 
led to the pioneering work on topological insulators in 
two and three spatial dimensions.^ ^ 

Another interesting gapped phase on the honeycomb 
lattice is the Kekule phase, illustrated in Fig. lb. The 
Kekule phase is topologically trivial but it is character- 
ized by a Z3 order parameter. The latter describes three 
degenerate Kekule ground states, obtained by translat- 



<l c 

1 

1 




1 9 

3 


' 4^ 


K — 7\ 
C — ^ 




4 


1 c 


■2 y 

) <• 




(a) 



(b) 



FIG. 1: The model: (a) Square lattice with |$o magnetic 
flux per plaquette and dimerized hopping amplitudes. The 
± on the left for each row show the choice of gauge for the 
Peierls phase factors. The thick (thin) bonds indicate an in- 
creased (decreased) hopping amplitude in the x (green) and 
y (blue) directions, the •/o dots on the lattice sites represent 
an increase/decrease in local charge density, and the dotted 
red line indicates the NNN hopping (showing only inside a 
single plaquette). The four sites of the unit cell are marked 
1-4. (b) Honeycomb lattice with similar color scheme 



ing the pattern depicted in Fig. lb by the Bravais lat- 
tice primitive vectors. Singly quantized vortices in the 
phase of this order parameter have been shown to possess 
a stable fermionic zero mode, carry fraction al cha rgj^ 
±e/2, and obey fractional exchange statistics .^i^Ell Simi- 
lar physics is realized in the dimerized phase on the 7r-flux 
square lattice.'^ 

In order to map out these interesting phases, Raghu 
et aZ.'i^ studied a simple model describing fermions, both 
spinless and spinful, on the honeycomb lattice with short 
range interactions. They found phase diagrams contain- 
ing CDW, SDW, QAH and QSH order but did not con- 
sider the possibility of the Kekule phase. With the goal of 
complementing this previous work we study similar mod- 
els and find that the Kekule phase is in fact present over 
a large portion of the phase diagram. Specifically, for 
spinless fermions, Kekule order is established when the 
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nearest-neighbor (NN) and next-nearest-neighbor (NNN) 
repulsion parameters Vi and V2 are of similar size. On the 
TT-flux square lattice we obtain a phase diagram with the 
Kekule phase replaced by a simpler dimerized phase hav- 
ing similar properties. In addition, we find stripe phases 
that tend to dominate over the QAH and the dimerized 
phases when V2 is large. Nevertheless, when the third- 
neighbor repulsion V3 is included, the Kekule and QAH 
phases are stabilized over a portion of the phase dia- 
gram. For spinful fermions, the possibility of a novel 
spin-Kekule phase characterized by the emergence of spa- 
tially modulated spin-dependent hopping amplitude is 
explored. 

Having a mean-field theory of the dimerized phase for 
the TT-fiux model allows us to perform a self-consistent 
calculation of the electronic structure for a system hav- 
ing a vortex in the dimer order parameter, which is set 
up in a fashion similar to Ref. 1171 The underlying Z4 
symmetry of the lattice raises the issue of whether it is 
reahstic to talk of J7(l) vortices in such a system, but as 
we have shown previously the 4-fold anisotropy induced 
by the lattice is very weak and should support such vor- 
tices on length scales long compared to the lattice spac- 
ing. Further evidence is supplied here to support this 
claim. 

II. SPINLESS FERMIONS 

We begin with the tight-binding Hamiltonian for spin- 
less fermions with NN and NNN interactions, H = 
Ho + Hj, with 

Ho = -tJ2ie'''^ 4c, +h.c.), (1) 

and 

Hj^Vi^ n.rij +V2^ rurij. (2) 

Here cj is the creation operator for a fermion on site of 
the square or honeycomb lattice, the Peierls phase factor 

27r 

..,.-y^A..i (3) 

is defined on a link (ij), t is the hopping amplitude be- 
tween NN sites, Vi and V2 are the energy scales for NN 
and NNN interactions respectively and rii = c\ci is the 
number operator. We focus on Hamiltonians Hq that 
produce Dirac-type spectra for fermions in the absence 
of interactions. 

In the following we treat Hj in the mean- field (MF) 
approximation which should provide reliable information 
about the possible gapped phases in the phase diagram. 
We consider both the on-site and bond MF decoupling 
channels, 

niUj nt{nj) + nj{ni) - {jii){nj) (4) 
n^Tij -A.jclcj ~ A*jclc, + AijA*^ (5) 



where = (cjci) with i,j belonging to NN and NNN 
bonds. 



A. TT-fiux square lattice 

The system has a half magnetic flux quantum $0 = 
hc/e per plaquette. In this case the assumption of spin- 
less fermions is quite natural because within the context 
of the realization described in Ref. 2 the electron spins 
would be polarized along the field. As in Ref. [171 we 
choose the Landau gauge A — (<i>o/2)(— y, 0) where we 
have set the lattice spacing to unity and work at half fill- 
ing. The unperturbed Hamiltonian then has a spectrum 



E^' ^±2ty sin^ -I- cos2 fcy , (6) 

with Dirac points at (0,±7r/2). 

The repulsive interactions can produce semi-metal 
(SM), charge density wave (CDW), quantum anomalous 
Hall (QAH), stripe and dimerized phases. The unit cell 
for this set up has four basis sites and is displayed in Fig. 
[iji. The SM phase corresponds to the undistorted 7r-flux 
lattice, which has no gap in the spectrum. 

The CDW corresponds to a modulation in charge den- 
sity in both the x and y directions resulting in a checker- 
board pattern on the lattice, whereas the stripe phase 
will be modulated only along one direction, hence giving 
rise to 'stripes' of increased/decreased charge density. We 
use the following ansatz, 

{n,) = ^- + p{-iy^+^y + u{-iy^ (7) 

for these two phases where = {ixiiy)- 

The QAH phase, which is characterized by a quantized 
Hall conductance without Landau levels, occurs when a 
gap opens due to the spontaneous breaking of time rever- 
sal invariance. We require here that its order parameter, 
Aij = Ax+y, produces a NNN hopping consistent with a 
half-flux quantum per unit cell. 

The dimerized phase, with its order parameter Aij — 
Ax^y along NN this time, will act to increase/decrease the 
existing hopping amplitude t. It can nominally take on 
any value locally, but we concentrate on either constant 
dimerization or one that will support the quasi U(l) vor- 
tices mentioned above. We also note the energy scale t 
becomes isotropically renormalized for any flnite value of 
Vi yielding a contribution 5t. 

Introducing the Fourier transform, Cj = 
Ck, our Hamiltonian can then be 
brought into the matrix form in momentum space 

H = Y. ^k^k^-k + Eo (8) 

k 
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where = (c^^, 4k, 4 

Eo = N[ 



3k' C4k ) , and 



P 



8(^1-^2 -14s) 8(^2-1^3) 

e 



Vi 



^ ^2 , 2_^^2n 



(9) 



1^2 ' 14 

with p = 4p(V^2 + 1^3 - Vi), V ^ Av{V3 - V2), Vx = ViAi, 
r]y — ViAy and ^ = Ax:+yV2. The Hamiltonian matrix 
reads 



with 



Hk = 



r* p + D —iiy — ri* 



-a* 



r 



V p-" 
\-^y -^x r* -p-vj 



2(tcos fcajjiy + i?7a;,i/ sin kx^y), 



(10) 



(11) 



r = 2i^[cos(fc2, + fcj^) — cos{ky — kx)] (12) 

and t = t + 6t. For reasons that will become apparent 
shortly we have also included in the above analysis the 
third-nearest-neighbor (NNNN) repulsion V3. 

We can now diagonalize Hk and obtain the exact dis- 
persion (setting the chemical potential to zero) with four 
particle-hole symmetric branches 

+ IQ^^s^s^ ± 2(p2p2 + AD^Pcl - 32Dif]x£,Cyslsy (13) 

I A-2 2 2 I ie-2,'2 2 2 , e ,( (^2 2 2/ 2 2, 2 2\\3l5 

-H4i/ ?7yS,y-f 16p ^ s^Sy+6^ s^Sy{ri^s^ + r)ySy)}''\\ 

Here k is taken over the reduced Brillouin zone — f < 
kx ^"2 7 2^ — — 2^ and c^; ^ cos kx^ Sy ^ sin ky , etc. 
From here, we may calculate the free energy 



So 



(14) 



where /3 = 1/ feT and the summation is over all four 
branches of i?k- The ground state is then obtained by 
minimizing F with respect to each of the order param- 
eters. This yields a set of gap equations (listed in Ap- 
pendix A) from which the phase diagram (at T = 0) seen 
in Fig. [2^ below follows after a numerical self-consistent 
iteration. The transition from the SM phase to the CDW 
is second order, whereas all other transitions are first or- 
der. 

We observe that, within the present model with Vi and 
V2, the stripe instability completely wipes out the QAH 
and the dimerized phases expected to occur at finite V2. 
In order to suppress the stripe phase we include a fur- 
ther term in the Hamiltonian above, namely a V3 term 
for NNNN interactions. Such term will generically be 
present in any system that can support NN and NNN 
interactions. To simplify matters we shall ignore any 
contribution to the kinetic energy generated by V3 and 
focus on the frustration it generates for the stripe insta- 
bility. We find that for V3 non-zero both the QAH and 
the dimerized phases are stabilized (Fig. ^p) although 
stripes reappear at stronger V2. 
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FIG. 2: Phase diagram for 7r-flux model with spinless 
fermions. (a) Vi = 0. (b) V3 = 1.5t 



B. Honeycomb lattice 

On the honeycomb lattice we may set Oij = (no mag- 
netic field needed for Dirac band structure). The unper- 
turbed spectrum reads 



4°) = ±i|r(k)| 



(15) 



where r (k) = g*'^ '*' , with vectors shown in Fig. 
[TJd. Also, there is no option of including a stripe phase, 
leaving us with the CDW, Kekulc and QAH phases. As 
seen in Fig. [ija the system now has 6 atoms in the unit 
cell. 

For the CDW we set 



» 



(16) 



where A and B refer to the two sublattices of the undis- 
torted honeycomb lattice, which leads to the contribution 
(each site having 3 NN and 6 NNN sites) 



JiJcDW = 3p (I/i - 2^2) ^ {-iy + ^Np" {Vi - 2V2) 

. 

For the Kekule phase, we again use the decoupling in 
equation (jZj) where c[cj{(^-Ci) simply adds to the tc[c.j 
term in Hq according to 

t^ta = t + 5t + Viria, a =1,2, 3 (18) 
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and 



N 



(19) 



The index a above labels the three bonds emanating from 
each A site. We may also write 



2tt 

rja = rj COS [ip+ —a 



(20) 



where tp parametrizes the most general Kekule distortion. 
We shall find below that the ground state energy is min- 
imized by ip — 0, 27r/3, 47r/3, corresponding to the three 
degenerate Kekule ground states. For any (p it holds that 
J2Va^0, and J^vl 



For the QAH phase we have (cjcj) = — ±i£,a (for 
i,j G NNN) where a = A,B, and the sign is taken ac- 
cording to the arrow in Fig. [ij). This generates a contri- 
bution 

SHqah = ^2 E (^^^^I^J + ^-c-) + 3^^2(ei + ^1) (21) 

m) 

to the mean-field Hamiltonian. 

Combining these contributions and Fourier transform- 
ing we arrive at the fc-spacc Hamiltonian of the form Eq. 
(|8| with ^'k a six-component spinor and 
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Here ia = P = 3p {Vi - 2V2), h = ■ k, 



Sk = -i 



i(k3 



and 



Eq 



N 



P T_ 
3(Vi-2V2) Vi 



i{ki 



,6' + fB' 



(22) 



(23) 



V2 



St' 



. (24) 

The reduced Brillouin zone, which is derived below in the 
Appendix B, is comprised of the following set of points 



(fci,fc2,/c3) 



1,2, 



_ 2n 1 
L and 



(71 — m, m, —n) 



(25) 



iV is the total number 



where m,n 
of sites. 

It is possible to find exact eigenvalues of the 6x6 ma- 



trix indicated in Eq. (|22| by noticing that is block- 



diagonal with two 3x3 blocks on the diagonal. Unfor- 
tunately, the eigenvalues obtained as roots of the asso- 
ciated cubic secular equations are given by lengthy and 
complicated expressions, which do not lend themselves 
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FIG. 3: Phase diagram for honeycomb lattice. At the mean- 
field level all transitions from the SM phase are second order 
whereas transitions between all the gapped phases are first 
order. The x's represent the line along which the ratio V2/V1 
would be expected to fall in graphene based on a crude esti- 
mate of the bare Coulomb repulsion.'^ 



to a convenient analysis. For this reason we choose to 
calculate the eigenvalues of TYk numerically for a dense 
discrete set of momenta k in the first BZ. We use these 



eigenvalues to compute the free energy, Eq. ( 14 1 , for a set 
of mean field parameters p, fj, (p and and given fixed 
values of Vi and V2 ■ The ground state of the system is 
then found by minimizing the free energy (at T = 0) with 
respect to the above MF parameters. This is achieved by 
employing the standard Powel multivariate minimization 
routine. 

The resulting phase diagram is displayed in Fig.[3j We 
observe that in addition to the CDW and QAH phases 
identified previously^* a large portion of the phase dia- 
gram is occupied by the Kekule phase. 

Although the critical lines in the phase diagram are de- 
termined numerically, the critical couplings for the CDW 
at V2 = and for QAH at Vi = can be found exactly 
from the corresponding gap equations for the model ex- 
cluding Kekule order. These read 



t 

Wc 

t 



3 

N 



y— 

^ \t(\ 



(k) 



1.341, 



2 v-^ (yi sink • b,) 

E , s , =i 0.840, 



3iV 



|r(k) 



(26) 
(27) 



where the correspond to the set of NNN vectors for 
cither of the triangular sublattices inside a single plaque- 
ttc. If we wish to express critical couplings in terms of 
the bare hopping amplitude i we must include the MF 
equation for the hopping renormalization 



1 



'^^ = ^i^E 1^(^)1^ 0-262^1 



(28) 
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We thus obtain Vic ^ 0.93^ and V2c ^ l-20t, in good 
agreement with the numerical results of Fig. [3] Our value 
of also agrees with Ref. [TSJ but our Vic is smaller 
by a factor of about 1.5 when expressed in terms of t 
and about a factor of 2 in terms t. We do not know 
what is the reason for this discrepancy. Since the critical 
couplings quoted above agree with our numerically de- 
termined phase diagram, we are confident that these are 
correct within the definitions employed in this study. 



III. SELF-CONSISTENT VORTEX STRUCTURE 

In a superconductor or a superfluid the order param- 
eter has a global U(l) symmetry related to its complex 
phase. This means that, even on a lattice, a U(l) vortex 
is a legitimate topological defect. In the Kekule (dimer- 
ized) phase the order parameter exhibits a global Z3 (Z4) 
symmetry leading to the possibility of a Z3 (Z4) vortex, 
which can be pictured as a point where the correspond- 
ing 3 (4) domains meet. Since the domain walls cost an 
energy per unit length, the energy of a single isolated Z„ 
vortex diverges linearly with the system size and, equiv- 
alently, a vortex-antivortex pair experiences linear con- 
finement. U(l) vortices, on the other hand, experience 
only logarithmic confinement. This has important conse- 
quences for the appearance of vortices in these systems. 
U(l) vortices can thermally unbind above the Kosterlitz- 
Thouless transition temperature T^t while Z„ vortices 
remain confined at all temperatures. 

It has been suggested^^ that the Z„ vortices in the 
Kekule/dimerized phase can nevertheless be observed 
since at short length scales they resemble U(l) vortices, 
and thus, for relatively short intervortex separations, in- 
teract only logarithmically. This has to do with the fact 
that the energy cost of a domain wall is relatively small. 
At length scales exceeding the confinement length Cconf, 
vortices remain linearly confined, as dictated by symme- 
try. If Cconf is sufficiently long compared to the vortex 
core size (^y, however, then the zero mode and the frac- 
tional charge associated with an individual vortex could 
be observed experimentally. To address this issue quan- 
titatively we now carry out a fully self-consistent calcula- 
tion of a vortex in the dimerized phase on the 7r-ffux lat- 
tice. Our calculations show that, indeed, at short length 
scales, a vortex in the dimer order parameter resembles a 
U(l) vortex while on longer length scales clear domains 
separated by domain walls emerge. 

Within our region of parameter space on the square 
lattice having a stable dimerized phase, we set up a dis- 
cretized version of a U(l) vortex in the dimerization pat- 
tern Aij as described in Ref. IT. For this initial vortex 
the corresponding MF Hamiltonian 



using the condition 

A,, = {4c,) = UlU.kidjdk) - Yl ^lU,if{ei) (30) 

I 



where / is the fermi function, and e; and dj are the eigen- 
values and the eigenmodes of Hmf- The unitary matrices 
Uij are comprised of the eigenvectors of the correspond- 
ing Hamiltonian matrix Ti-ij, which we find numerically 
through exact diagonalization. Eqs. ( 29|30 1 are then it- 
erated to self-consistency. We have solved the system 
for odd lattice sizes up to 49 x 49 using open boundary 
conditions with a vortex positioned at the central site. 

Our results for the self-consistent vortex structure are 
summarized in Fig. li] Instead of (real) bond fields Aij 
it is useful to consider a complex on-site 'dimer' order 
parameter defined as 



9i 



llA+nA. 



i.i-\-l.i^i,i-\-fl- 



(31) 



Here the jl are the four nearest-neighbor unit vectors, 
li,i+x = (-1)'", Ji,i+y = i(-l)''' and 7i,i-/i = 7i-/i,i- 
The phase of gi contains information on the local ori- 
entation of the dimer pattern, e.g. gt ~ e*5^"^2) with 
n = 1,2,3,4 describes four basic 'columnar' patterns 
modulated along the x and y directions while gi ^ e*^^" 
describes four possible 'box' phases with the strength of 
dimerization modulated along both the x and y direc- 
tions. We note that in the uniform system the box phases 
correspond to the ground states. 

The vortex has a finite size core with a radius inversely 
proportional to the gap A, as seen in Fig. |4Jd. Near the 
core the phase behaves as in a U(l) vortex although at 
longer distances domain walls are seen to form along the 
lattice diagonals, Fig.|4^, characteristic of a Z4 vortex. A 
net fractional charge e/2 accumulates near the co re. Fig . 
|4|:, as expected on the basis of general arguments .^3111 

The nature of the dimerized box ground state and the 
competition with other possible states can be understood 
from the following analysis. If we tune Vi and V2 so 
that we are in the dimerized phase, then all other or- 
der parameters can be set to zero. In the uniform state 
we can parameterize rjx and rjy by an angle x, so that 
A — {rjx, rjy) = A(cos X, sin x), and write our ground 
state energy per site as (summing over all negative en- 
ergy states) 



Eoix) = -E' 



(32) 



with 



A2 



= 4t (cos + cos ky) + — (sin k^ + sin ky), 



Wj. = cos 2x(sin fc^, — sin ky). 



Hmf = Ho + YiA,jclcj + h.c.) (29) Taylor expanding 



is diagonalized and the new order parameter is found 




(33) 
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FIG. 4: Self-consistent vortex structure in the dimerized 
phase, (a) Phase and (b) magnitude of the on-site dimer 
order parameter Qi defined in Eq. (31 1, (c) Charge density 



profile showing excess charge density in the vortex core. The 
total accumulated charge is e/2. (d) Anisotropy parameter 6 
defined in Eq. (|35|. For figures (a)-(c), A/t ^ 0.24 



(sin^ krc — sin^ ky)^ 



wc then have to fourth order in tUk 



k 



32 



cos' 



2xE 



(34) 

The ground state is clearly minimized by x = for 
n = 1,3, corresponding to the box phase. The columnar 
phase corresponds to the maximum of energy. However, 
by plotting the anisotropy ratio 



6 = 



EojO) - Eo{Tr/A) 
Eoin/A) 



(35) 



against A/t, one can see, as in Fig.|4|i, that the difference 
between the maximum and minimum values of i?o(x) "^iU 
be very small for a broad range of values of A/t. The 
smallness of 6 explains why at short length scales our Z4 
vortex behaves as a U(l) vortex. 



IV. SPINFUL FERMIONS 

For spinful Fermions a Hubbard term C/ '^^T'^U 
must be added to the interaction Hamiltonian, Eq. ([2|, 
reflecting the strong on-site Coulomb repulsion between 
electrons of opposite spin. It is well-known that for a 
bi-partite lattice, a large U favors the antiferromagnetic 
SDW state. The mean-field phase diagram for spinful 
electrons on the honeycomb lattice in the U-V1-V2 space 
has been mapped out in Ref. l18] For small values of U 



the phase diagram resembles that shown in Fig. [3]except 
that the QAH phase is replaced by the QSH phase. The 
two are degenerate at the mean field level but quantum 
fluctuations favor the latter. 

A question that we would like to answer is how the 
U-V1-V2 phase diagram on the honeycomb lattice is af- 
fected by the presence of the Kekule phase, which was 
not considered in Ref. 1181 Is there a region in the param- 
eter space where a 'spin Kekule' phase is stabilized? The 
latter represents a spinful version of the ordinary Kekule 
phase and its order parameter is deflned as follows. The 
nearest neighbor interactions in Hj can be re- written us- 
ing the identity {n, - l){nj - 1) = 1 - i (xf^)^ xfj, where 
Xij = 4a'^ai3^jl3' M = • ■ • 3 and CT^ = (1, cr). In a mean 
field picture, then, (x°) ^ corresponds to the standard 
Kekule phase and (x*) 7^ describes the spinful version. 

We find, however, that the spin Kekule phase is not 
a mean field ground state of our Hamiltonian. Instead, 
a state in which two projections of the spin form the 
ordinary Kekule phase is favored. This can be seen from 
the following symmetry based argument. We focus on 
the /i = 3 state and put the spin up into the ordinary 
Kekule phase with parameters (77, ip) that minimize its 
ground state energy. The spin down electrons must 
then go into the Kekule phase with parameters (— ??, f) 
due to the (T3 in the order parameter. However, this 
is not a ground state for spin down electrons, meaning 
that the overall energy must be higher than the /Lt = 
state. Thus the energy of the SK phase is not exactly 
degenerate with the standard Kekule phase as is the case 
for the QAH/QSH spectrum, but instead has a small 
splitting with the Kekule phase being favored as the 
ground state. So, in absence of some type of interaction 
that could favor the SK phases, our current Hamiltonian 
is not capable of producing this potentially interesting 
phase. 



V. CONCLUSIONS 

Our studies show that interacting fermions moving on 
the honeycomb and 7r-fiux square lattices can form a sta- 
ble Kekule phase, in addition to the previously identi- 
fied CDW, SDW, and QAH/QSH phases. Two questions 
naturally arise: Can the Kekule phase be stabilized in 
a realistic system and if so, how can it be experimen- 
tally observed? Inspection of the phase diagram shown in 
Fig. 2 reveals that the prospects of observing the Kekule 
(dimerized) phase on the square lattice are not very good. 
When the NNN interaction is sufficiently strong to sup- 
press the CDW, one obtains the stripe phase. Even 
when this is suppressed by introduction of NNNN re- 
pulsion, the dimerized phase appears only at strong cou- 
pling (Vi ~ 6t) and only in a small sliver of the parameter 
space. We thus conclude that although the square lattice 
is very convenient from the point of view of theoretical 



7 



considerations, it is an unlikely candidate for the experi- 
mental observation of the Kekule/dimerized phase. 

Prospects for the Kekule phase appear much better 
on the honeycomb lattice as it occupies a large portion 
of the V1-V2 mean field phase diagram shown in Fig. 3. 
The Kekule phase should appear whenever the interac- 
tion strength becomes large and Vi, V2 are comparable. 
In natural graphene the Coulomb interaction is strong 
but evidently not strong enough to open up a gap at 
the Dirac point P at least when a sample is placed on a 
substrate such as Si02. It has been suggested recently, 
on the basis of large scale Monte Carlo simulations, that 
graphene freely suspended in vacuum could become an 
insulator when the Coulomb interaction between elec- 
trons is accounted for.^^ The difference between the sub- 
strate and vacuum arises from the dielectric constant, 
which weakens the Coulomb interaction in the vicinity 
of a polarizable solid such as Si02, and also because of 
disorder induced by the substrate. Although Ref. .19. fo- 
cused on the CDW instabihty, one may argue based on 
our present mean field study that the Kekule phase is a 
more likely candidate for an instability induced by the 
Coulomb interaction. In a very crude estimate the ratio 
V2/V1 should be equal to the ratio between the NN and 
NNN bond lengths, i.e. V2/V1 w l/\/3 ~ 0.5773. As il- 
lustrated in Fig. 3, this ratio of V2/V1 clearly favors the 
Kekule phase over CDW if the overall effective interac- 
tion strength is strong enough to open up a gap. Recent 
experimental observation of the fractional quantum Hall 
effect in suspended graphene^*^ confirms the increased im- 
portance of interactions over graphene on a substrate, 
where the fractional quantum Hall effect thus far eluded 
experimental detection. At zero field, however, even the 
suspended samples appear semi-metallicP^ down to 1.2K, 
suggesting that Coulomb interactions are still too weak 



to open up a significant gap. We note, finally, that the 
Kekule phase has been proposed to emerge as the lead- 
ing instability of graphene in an applied magnetic field.l^il 
Also, a phase analogous to the Kekule phase in graphene 
has been proposed to exist for electrons on the kagome 
lattice at 1/3 filling.l22l 

The Kekule phase has broken translational symmetry 
characterized by a wavevector connecting the two in- 
equivalent Dirac points in the graphene spectrum. The 
on-site charge density remains uniform and this limits 
the spectrum of experimental probes capable of detecting 
this pattern of symmetry breaking. Here we propose to 
use the technique of Fourier transform scanning tunnel- 
ing spectroscopy (FT-STS) that has already been applied 
to graphene. ^■^ With sufficient resolution, FT-STS allows 
the mapping out of fine details of quasiparticle interfer- 
ence patterns at non-zero momenta, and, with help from 
the theoretical modeling of such patterns,'^ it should be 
possible to establish the existence of the Kekule or other 
symmetry breaking phases. 
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VI. APPENDIX 

A. Gap Equations 

For completeness we list the gap equations for the tt- 
flux model on the square lattice in the limit T ^ for P, 
P, Vx, Vy I C and 5t : 



P = 



Vy 



St- 



' --^^ E ^ ± ^ i^'P + ^^Pe K sin^ ky)] 



4 (V2 V3) ^ — l^j^ J_ ^^p2 ^^p ^^g2 _ ]^g^^^ COS ky sin^ kx sin ky + 4.vrjy sin^ ky)'] 



2Vi 



AT E — i" ['^^ '^^^ kx± -t^ sw? kx (l6^^rya; sin^ k^ sit? ky - Avi£^ cos ky sin ky)] 
fe,± Ek 



21/1 



~ E TrTTT ^y^^ '^^^^y ^"^^^ ^ '^^(^'ny sin"* ky sw? kx)] 



i = E [£, sin^ ky±^ [p^i s\T? ky + 4^(77^ sin^ kx sit? ky + rj^ sin'' ky) - vtrjx cos ky sin ky)] 

I. J- \Ek.± ^■ 



fe,± 



Ek 



r^r^ — r P (cos^ ky -|- cos^ k^ ± ^ (2i^^i cos^ ky — Si'i'qx^ cos ky sir? kx sin ky)] 
k,± I '^'il 



(36) 
(37) 
(38) 
(39) 
(40) 
(41) 
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where 

Ek — (p'^P^ + ^v^'P cos^ ky — 32pfr/a;f cos ky sin^ k^ sin ky + Av'^rj^ sit? ky 

+ l&p^i^ SIT? k^ SIT? ky + 64^^ sin^ k^, sit? ky{-ql sit? k^ + rj^ sit? ky)) ' (42) 



As mentioned above, the phase diagrams seen in Fig. 
[2] were mapped out by solving these equations self- 
consistently over the desired region of parameter space. 
That is to say, an initial set of mean field values were 
substituted into the equations and then checked against 
the LHS of each equation. If the difference was greater 
than some tolerance set at the beginning, the new values 
were fed back in until convergence was achieved. 



B. Brillouin Zone 

The underlying bravais lattice is spanned by the 
primitive vectors (3ai,3a2) = (Ai,A2) where Ai = 
I (\/3i + y) 0,0, A2 = 3yao and aq is the lattice spac- 



ing. The reciprocal lattice vectors, Gi 



27r 2 ' 



and 



G2 — 1^ {V^y ~ 2;) , then follow and we have for the 
brillouin zone 



mGi 
27r 



- nG2 
1 r 



3y/3ao L 



Thus with ai = 

a; + I y j flQ and as 
/c2 = a2 • k 



TO, n e N 
(2to — n) X + ^/Sny 



V3a 



27T 71 — m 



and fc; 



(43) 



{^-^x+ ^yj ao, a2 

yao we get, fci = ai • k 
as • k 



2tt_ n 
3 L- 
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